!dk budget
      subroutine budget 
!-----------------------------------------------
!   M o d u l e s 
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double 
      use corgan_com_M 
      use cindex_com_M 
      use cophys_com_M 
      use blcom_com_M
      use numpar_com_M 
      use objects_com_M, ONLY: link_obj, iphead_obj 
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02  
!...Switches: -yf -x1             
!                                                       *
!********************************************************
!
      implicit none
!-----------------------------------------------
!   G l o b a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer , dimension(8) :: iv 
      integer , dimension(4) :: lioinput, lioprint 
      integer :: jsurf, n, nplost_neg, nplost_pos, np, is, newcell, newcell_nxt&
         , nptotl_obj 
      real(double), dimension(1) :: ixi1, ixi2, ieta1, ieta2 
      real(double) :: ixi4, ieta4, ixi5, ieta5, mupar, jxc, jyc, jtx, jty, tee&
         , toue, toui, tex, tey, tez, totalj, totarea_y, xc, yc, zc, phi, theta&
         , r_minor, tbex, tbey, tbez, totlchg, eb, qlost_neg, qlost_pos, &
         qavg_neg, qavg_pos, tparke, tparki, tparmx, tparmy, tparmz, tparkex, &
         tparkey, tparkez, tparkix, tparkiy, tparkiz
      logical :: expnd, nomore 
      character :: name*80 
!-----------------------------------------------
!
!
      tee = 0.0 
      toue = 0.0 
      toui = 0.0 
      tex = 0.0 
      tey = 0.0 
      tez = 0.0 
!
!     calculate the total toroidal current
!
      jsurf = 2 
      call area_bc (2, kbp1, 2, ibp1, jsurf, kwid, iwid, jwid, 1., x, y, z, &
         area_x, area_y, area_z) 
!
      totalj = 0.0 
      totarea_y = 0. 
!
      do k = 2, kbp1 
         do i = 2, ibp1 
!
            ijk = (k - 1)*kwid + (jsurf - 1)*jwid + (i - 1)*iwid + 1 
!
            jxc = 0.25*(pixx(ijk)+pixx(ijk+iwid)+pixx(ijk+iwid+jwid)+pixx(ijk+&
               jwid)) 
!
            jyc = 0.25*(pixy(ijk)+pixy(ijk+iwid)+pixy(ijk+iwid+jwid)+pixy(ijk+&
               jwid)) 
!
            xc = 0.25*(x(ijk)+x(ijk+iwid)+x(ijk+iwid+jwid)+x(ijk+jwid)) 
!
            yc = 0.25*(y(ijk)+y(ijk+iwid)+y(ijk+iwid+jwid)+y(ijk+jwid)) 
!
            zc = 0.25*(z(ijk)+z(ijk+iwid)+z(ijk+iwid+jwid)+z(ijk+jwid)) 
!
            if (rmaj > 0.) then 
               call toroidal_coords (pi, rmaj, xc, yc, zc, phi, theta, r_minor) 
            else 
               phi = 0. 
            endif 
!
!
            jtx = -jxc*sin(phi) 
            jty = jyc*cos(phi) 
!
            totalj = totalj + jtx*area_x(ijk) + jty*area_y(ijk) 
!
            totarea_y = totarea_y + area_y(ijk) 
         end do 
      end do 
!
!
      do n = 1, nvtx 
!
         ijk = ijkvtx(n) 
!
         toue = toue + 0.5*(jxs(ijk,2)**2+jys(ijk,2)**2+jzs(ijk,2)**2)*qdnv(ijk&
            ,2)/qom(2) 
         toui = toui + 0.5*(jxs(ijk,1)**2+jys(ijk,1)**2+jzs(ijk,1)**2)*qdnv(ijk&
            ,1)/qom(1) 
         tex = tex + ex(ijk)**2*vvol(ijk)*0.5 
         tey = tey + ey(ijk)**2*vvol(ijk)*0.5 
         tez = tez + ez(ijk)**2*vvol(ijk)*0.5 
!
      end do 
      tee = tex + tey + tez 
!
!
      tbe = 0. 
      tbex = 0. 
      tbey = 0. 
      tbez = 0. 
!
      totlchg = 0.0 
!
      do n = 1, ncells 
!
         ijk = ijkcell(n) 
!
         eb = 0.5*(bxn(ijk)**2+byn(ijk)**2+bzn(ijk)**2)*vol(ijk) 
         tbe = tbe + eb 
         tbex = tbex + 0.5*bxn(ijk)**2*vol(ijk) 
         tbey = tbey + 0.5*byn(ijk)**2*vol(ijk) 
         tbez = tbez + 0.5*bzn(ijk)**2*vol(ijk) 
!
         totlchg = totlchg + qdnc(ijk)*vol(ijk) 
!
      end do 
!
!
!     impose zero current conditions at boundary
!
      nplost_neg = 0 
      nplost_pos = 0 
      qlost_neg = 0.0 
      qlost_pos = 0.0 
!
      np = iplost 
!
  210 continue 
      if (np > 0) then 
!
         if (qpar(np) < 0.0) then 
!
            qlost_neg = qlost_neg + qpar(np) 
            nplost_neg = nplost_neg + 1 
         else 
            qlost_pos = qlost_pos + qpar(np) 
            nplost_pos = nplost_pos + 1 
         endif 
         np = link(np) 
         go to 210 
      endif 
!
      qavg_neg = qlost_neg/(float(nplost_neg) + 1.E-10) 
      qavg_pos = qlost_pos/(float(nplost_pos) + 1.E-10) 
!
      if (totlchg > 0.0) then 
         plost_neg = (plost_neg*(totlchg - qlost_pos))/(qlost_neg - 1.E-20) 
         plost_pos = 1. 
         plost_neg = amin1(plost_neg,1.) 
         plost_neg = amax1(plost_neg,0.) 
!
      else 
         plost_neg = 1. 
         plost_pos = (plost_pos*(totlchg - qlost_neg))/(qlost_pos + 1.E-20) 
         plost_pos = amin1(plost_pos,1.) 
         plost_pos = amax1(plost_pos,0.) 
!
      endif 
!
!
!     override ... set plost=0.
      plost_pos = 0. 
      plost_neg = 0. 
!
!     return discarded particles to the holding list
!
 2221 continue 
      if (iplost /= 0) then 
         np = iplost 
         iplost = link(np) 
         link(np) = iphead(1) 
         iphead(1) = np 
         go to 2221 
      endif 
!
!
!
 
      ijkctmp(:ncells) = ijkcell(:ncells) 
      iphd2(ijkcell(:ncells)) = 0 
!
!
      tparke = 0.0 
      tparki = 0.0 
      tparmx = 0.0 
      tparmy = 0.0 
      tparmz = 0.0 
      tparkex = 0.0 
      tparkey = 0.0 
      tparkez = 0.0 
      tparkix = 0.0 
      tparkiy = 0.0 
      tparkiz = 0.0 
!
      nptotl = 0 
      numtot(:nreg) = 0 
!
!
      newcell = ncells 
    1 continue 
      nomore = .TRUE. 
!
!dir$ ivdep
!
      newcell_nxt = 0 
      do n = 1, newcell 
         if (iphead(ijkctmp(n)) == 0) cycle  
         newcell_nxt = newcell_nxt + 1 
         ijkctmp(newcell_nxt) = ijkctmp(n) 
         nptotl = nptotl + 1 
      end do 
!
      newcell = newcell_nxt 
!
      if (newcell /= 0) then 
         nomore = .FALSE. 
!
         nplost = 0 
         np = iplost 
  601    continue 
         if (np > 0) then 
            nplost = nplost + 1 
            np = link(np) 
            go to 601 
         endif 
!
!
!     calculate the particle energy
!
         do n = 1, newcell 
!
            ijk = ijkctmp(n) 
            np = iphead(ijk) 
!
            if (np == 0) cycle  
!
            is = ico(np) 
!
!
!     electrons only (species 2)
!
!     if(is.eq.2) then
            if (is == 1) then
            tparki = tparki + qpar(np)*(up(np)**2+vp(np)**2+wp(np)&
               **2)/qom(is)
            tparkix = tparkix + qpar(np)*up(np)**2/qom(is) 
            tparkiy = tparkiy + qpar(np)*vp(np)**2/qom(is) 
            tparkiz = tparkiz + qpar(np)*wp(np)**2/qom(is)  
            end if
            if (is == 2) then
            tparke = tparke + qpar(np)*(up(np)**2+vp(np)**2+wp(np)&
               **2)/qom(is) 
            tparkex = tparkex + qpar(np)*up(np)**2/qom(is) 
            tparkey = tparkey + qpar(np)*vp(np)**2/qom(is) 
            tparkez = tparkez + qpar(np)*wp(np)**2/qom(is) 
            end if
            tparmx = tparmx + qpar(np)*up(np)/qom(is) 
            tparmy = tparmy + qpar(np)*vp(np)/qom(is) 
            tparmz = tparmz + qpar(np)*wp(np)/qom(is) 
!        end if
!
!     save total number in each region
            numtot(is) = numtot(is) + 1 
!
         end do 
!
         do n = 1, newcell 
            ijk = ijkctmp(n) 
            np = iphead(ijk) 
            if (np <= 0) cycle  
            iphead(ijk) = link(np) 
            link(np) = iphd2(ijk) 
            iphd2(ijk) = np 
!
         end do 
      endif 
!
      if (.not.nomore) go to 1 
!
      do n = 1, ncells 
!
         ijk = ijkcell(n) 
!
         iphead(ijk) = iphd2(ijk) 
!
         iphd2(ijk) = 0 
!
      end do 
!
!
!
 
      ijkctmp(:ncells) = ijkcell(:ncells) 
      iphd2(ijkcell(:ncells)) = 0 
!
!
!
      nptotl_obj = 0 
!
      newcell = ncells 
   11 continue 
      nomore = .TRUE. 
!
!dir$ ivdep
!
      newcell_nxt = 0 
      do n = 1, newcell 
         if (iphead_obj(ijkctmp(n)) == 0) cycle  
         newcell_nxt = newcell_nxt + 1 
         ijkctmp(newcell_nxt) = ijkctmp(n) 
         nptotl_obj = nptotl_obj + 1 
      end do 
!
      newcell = newcell_nxt 
!
      if (newcell /= 0) then 
         nomore = .FALSE. 
!
!
!
         do n = 1, newcell 
            ijk = ijkctmp(n) 
            np = iphead_obj(ijk) 
            if (np <= 0) cycle  
            iphead_obj(ijk) = link_obj(np) 
            link_obj(np) = iphd2(ijk) 
            iphd2(ijk) = np 
!
         end do 
      endif 
!
      if (.not.nomore) go to 11 
!
      do n = 1, ncells 
!
         ijk = ijkcell(n) 
!
         iphead_obj(ijk) = iphd2(ijk) 
!
         iphd2(ijk) = 0 
!
      end do 
!
!
!
      write (6, 1093) ncyc, nptotl_obj 
 1093 format(" ncyc=",i5," nptotl_obj=",i7) 
!
!
!     store history variables
!
!      tparke=tparke-tparmx**2-tparmy**2-tparmz**2
!      tparkex=tparkex-tparmx**2
!      tparkey=tparkey-tparmy**2
!      tparkez=tparkez-tparmz**2
      nh = mod(ncyc,nhst) 
      nh = max(1,min(nhst,nh)) 
      efnrg(nh) = tee 
      efnrgx(nh) = tex 
      efnrgy(nh) = tey 
      efnrgz(nh) = tez 
      ebnrg(nh) = tbe 
      ebnrgx(nh) = tbex 
      ebnrgy(nh) = tbey 
      ebnrgz(nh) = tbez 
      eoenrg(nh) = toue 
      eoinrg(nh) = toui 
      ekenrg(nh) = tparke*.5d0 
      ekenrgx(nh) = tparkex*.5d0
      ekenrgy(nh) = tparkey*.5d0 
      ekenrgz(nh) = tparkez*.5d0 
      ekinrg(nh) = tparki*.5d0
      ekinrgx(nh) = tparkix*.5d0
      ekinrgy(nh) = tparkiy*.5d0 
      ekinrgz(nh) = tparkiz*.5d0  
      charge(nh) = totlchg 
      toroidal_j(nh) = totalj 
      thistry(nh) = t 
!
!
      write (6, 1091) ncyc, nptotl, totlchg 
 1091 format(" ncyc=",i5," nptotl=",i7," totlchg=",1p,e12.5) 
      write (*, *) 'nplost', nplost 
!
!
      write (6, 1092) tparke + tparki, tee, tbe 
 1092 format(" total particle energy=",1p,e10.3," electric field nrg=",1p,e10.3&
         ," magnetic field nrg=",1p,e10.3) 
!      write (6, 1094) tparkex, tparkey, tparkez 
!      write (56, *) tparkex, tparkey, tparkez 
 1094 format(" total particle energy x=",1p,e10.3," y=",1p,e10.3," z=",1p,e10.3&
         ) 
      write (6, 1095) tex, tey, tez 
      write (57, *) tex, tey, tez 
 1095 format(" total field x=",1p,e10.3," y=",1p,e10.3," z=",1p,e10.3) 
      write (6, 1096) tparmx, tparmy, tparmz 
      write (58, *) tparmx, tparmy, tparmz 
 1096 format(" total momentum x=",1p,e10.3," y=",1p,e10.3," z=",1p,e10.3) 
!
      chcons1(nh) = chnorm1 
      chcons2(nh) = chnorm2 
!
      return  
      end subroutine budget 
